# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/7/24

library(ggrepel)
library(optparse)
library(magrittr)
library(ggpubr)
library(tidyverse)
library(FSA)
library(RJSONIO)
library(car)

createWhenNoExist <- function(f) {
  !dir.exists(f) && dir.create(f)
}

self.my.fc <- function(data1, data2, method) {
  fc <- 0
  curMethod <- method
  if (curMethod == "anova") {
    curMethod <- "mean"
  }else {
    curMethod <- "median"
  }

  if (curMethod == "median") {
    median1 <- median(data1)
    median2 <- median(data2)
    fc <- median2 / median1
  }else if (curMethod == "mean") {
    mean1 <- mean(data1)
    mean2 <- mean(data2)
    fc <- mean2 / mean1
  }
  if (is.nan(fc)) {
    fc <- 1
  }
  list(fc = fc, method = curMethod)
}

option_list <- list(
  make_option("--i", default = "AllMet_Raw.txt", type = "character", help = "metabolite data file"),
  make_option("--g", default = "group.txt", type = "character", help = "sample group file"),
  make_option("--sc", default = "sample_color.txt", type = "character", help = "sample color file"),
  make_option("--config", default = "calculate_config.json", type = "character", help = "config file")
)
opt <- parse_args(OptionParser(option_list = option_list))

args <- commandArgs(trailingOnly = F)
scriptPath <- dirname(sub("--file=", "", args[grep("--file", args)]))
source(str_c(scriptPath, "/base.R"))

configJson <- fromJSON(opt$config)

diffMethod <- configJson$method
pValue <- configJson$pValue %>%
  as.numeric()
fdr <- configJson$fdr %>%
  as.numeric()
fcMethod <- configJson$fcMethod
isInter <- configJson$isInter %>%
  as.logical()
execPostHoc <- configJson$execPostHoc %>%
  as.logical()
execTrend <- configJson$execTrend %>%
  as.logical()
mcInfoList <- configJson$mcInfos
groupSort <- configJson$groups
trendMethod <- configJson$trendMethod

my.test <- function(data, method, metaboliteTb) {
  p <- 0
  curMethod <- method
  outTd <- metaboliteTb
  if (curMethod == "auto") {
    sTest <- data %>%
      group_by(ClassNote) %>%
      summarise(p = {
        var <- var(value)
        if (var == 0) {
          0
        }else shapiro.test(value)$p.value
      }, .groups = "keep")
    bTest <- leveneTest(value ~ ClassNote, data)
    bp <- bTest$`Pr(>F)`[1]
    if (all(sTest$p > 0.05) && bp > 0.05) {
      curMethod <- "anova"
    }else {
      curMethod <- "kw"
    }
  }

  if (curMethod == "kw") {
    test <- kruskal.test(value ~ ClassNote, data = data)
    p <- test$p.value
    if (is.na(p)) {
      p <- 1
    }
    curMethod <- "kruskal.test"
    if (execPostHoc) {
      dt <- dunnTest(value ~ ClassNote, data = data)
      outTb <- dt %>%
        .$res %>%
        as_tibble() %>%
        rename(group = Comparison, `p adj` = P.unadj) %>%
        select(c("group", "p adj")) %>%
        gather("p adj", "Value", -group) %>%
        spread(group, "Value") %>%
        select(-c("p adj")) %>%
        rename_all(function(x) {
          x %>%
            map_chr(function(groupName) {
              groups <- groupName %>%
                str_split("-") %>%
                unlist() %>%
                str_trim()
              str_c("P.adj.", groups[1], ".", groups[2])
            })
        }) %>%
        rowwise() %>%
        do({
          result <- as_tibble(.)
          names <- colnames(result)
          tb <- names %>%
            map_dfc(function(name) {
              groups <- name %>%
                str_replace(., "P.adj.", "") %>%
                str_split("\\.") %>%
                unlist()
              group1 <- groups[1]
              group2 <- groups[2]
              group1Data <- data %>%
                filter(ClassNote == group1) %>%
                .$value
              group2Data <- data %>%
                filter(ClassNote == group2) %>%
                .$value
              fc <- self.my.fc(group1Data, group2Data, curMethod)$fc
              log2FC <- log(fc, 2)
              pColumnName <- str_c("P.adj", ".", group1, ".", group2)
              fcColumnName <- str_c("FC", ".", group1, ".", group2)
              log2FcColumnName <- str_c("log2FC", ".", group1, ".", group2)
              p <- result[1, pColumnName] %>%
                as.numeric()
              tibble(p = p, fc = fc, log2Fc = log2FC) %>%
                set_colnames(c(pColumnName, fcColumnName, log2FcColumnName))
            })
          tb
        }) %>%
        ungroup()
      outTd <- outTd %>%
        bind_cols(outTb)
    }
  }else if (curMethod == "anova") {
    test <- aov(value ~ ClassNote, data = data)
    p <- summary(test)[[1]][[1, "Pr(>F)"]]
    curMethod <- "anova"
    if (execPostHoc) {
      td <- TukeyHSD(test)
      outTb <- td %>%
        .$ClassNote %>%
        as.data.frame() %>%
        rownames_to_column("group") %>%
        as_tibble() %>%
        select(c("group", "p adj")) %>%
        gather("p adj", "Value", -group) %>%
        spread(group, "Value") %>%
        select(-c("p adj")) %>%
        rename_all(function(x) {
          x %>%
            map_chr(function(groupName) {
              groups <- groupName %>%
                str_split("-") %>%
                unlist()
              str_c("P.adj.", groups[2], ".", groups[1])
            })
        }) %>%
        rowwise() %>%
        do({
          result <- as_tibble(.)
          names <- colnames(result)
          tb <- names %>%
            map_dfc(function(name) {
              groups <- name %>%
                str_replace(., "P.adj.", "") %>%
                str_split("\\.") %>%
                unlist()
              group1 <- groups[1]
              group2 <- groups[2]
              group1Data <- data %>%
                filter(ClassNote == group1) %>%
                .$value
              group2Data <- data %>%
                filter(ClassNote == group2) %>%
                .$value
              fc <- self.my.fc(group1Data, group2Data, curMethod)$fc
              log2FC <- log(fc, 2)
              pColumnName <- str_c("P.adj", ".", group1, ".", group2)
              fcColumnName <- str_c("FC", ".", group1, ".", group2)
              log2FcColumnName <- str_c("log2FC", ".", group1, ".", group2)
              p <- result[1, pColumnName] %>%
                as.numeric()
              tibble(p = p, fc = fc, log2Fc = log2FC) %>%
                set_colnames(c(pColumnName, fcColumnName, log2FcColumnName))
            })
          tb
        }) %>%
        ungroup()
      outTd <- outTd %>%
        bind_cols(outTb)
    }
  }
  list(p = p, method = curMethod, outTd = outTd)
}

options(digits = 3)

sampleInfo <- read_tsv(opt$g) %>%
  rename(SampleID = Sample) %>%
  select(c("SampleID", "ClassNote")) %>%
  mutate(ClassNote = factor(ClassNote, levels = unique(ClassNote)))

sampleIds <- sampleInfo$SampleID

groups <- unique(sampleInfo$ClassNote)

allData <- read_tsv(opt$i) %>%
  rename(Metabolite = 1) %>%
  select(c("Metabolite", sampleIds)) %>%
  rowwise() %>%
  do({
    result <- as_tibble(.)
    metaboliteTb <- result %>%
      select(c("Metabolite"))
    kwData <- result %>%
      select(-c("Metabolite")) %>%
      gather("SampleID", "value") %>%
      inner_join(sampleInfo, by = c("SampleID"))
    for (group in groups) {
      sampleIds <- sampleInfo %>%
        filter(ClassNote == group) %>%
        .$SampleID
      groupData <- result[sampleIds] %>%
        unlist()
      result[, paste0(group, ".Mean")] <- mean(groupData)
      result[, paste0(group, ".Median")] <- median(groupData)
      result[, paste0(group, ".SD")] <- sd(groupData)
      iqr11 <- quantile(groupData, 0.25)
      iqr12 <- quantile(groupData, 0.75)
      iqr1 <- paste0("[", iqr11, ",", iqr12, "]")
      result[, paste0(group, ".IQR")] <- iqr1
    }
    rs <- my.test(kwData, diffMethod, metaboliteTb)
    result$P <- rs$p
    mc <- if (rs$method == "anova") {
      "tukey"
    }else "dunn"
    result$test.method <- if (execPostHoc) {
      str_c(rs$method, ".", mc)
    }else {
      rs$method
    }
    result %>%
      left_join(rs$outTd, by = c("Metabolite"))
  }) %>%
  select(-c(sampleInfo$SampleID))

data <- allData %>%
  mutate(FDR = p.adjust(P, method = "fdr")) %>%
  select(-c("test.method"), everything()) %>%
  select(-c(starts_with("P."), starts_with("FC."), starts_with("log2FC.")), everything())

data

parent <- paste0("./")
allFileName <- paste0(parent, "/AllMet_Test.csv")
write.csv(data, allFileName, row.names = FALSE)

finalDf <- read_csv("AllMet_Test.csv") %>%
  rowwise() %>%
  do({
    result <- as.data.frame(.)
    p <- result[1, "P"] %>%
      as.numeric()
    q <- result[1, "FDR"] %>%
      as.numeric()
    value <- if (p < pValue & q < fdr) {
      1
    }else 0
    result$IS_Uni_P_Sig <- value
    result
  }) %>%
  ungroup() %>%
  filter(IS_Uni_P_Sig == 1) %>%
  select(-c("IS_Uni_P_Sig"))

finalDf <- if (execPostHoc) {
  finalDf %>%
    rowwise() %>%
    do({
      result <- as_tibble(.)
      bs <- mcInfoList %>%
        rowwise() %>%
        do({
          inResult <- as_tibble(.)
          groups <- inResult["group"] %>%
            str_split("-") %>%
            unlist()
          p <- inResult["p.adj"] %>%
            as.numeric()
          log2FC <- inResult["log2Fc"] %>%
            as.numeric()
          pColumn <- str_c("P.adj.", groups[1], ".", groups[2])
          log2FcColumn <- str_c("log2FC.", groups[1], ".", groups[2])
          pV <- result[1, pColumn]
          logFc <- result[1, log2FcColumn]
          b <- pV < p && !is.na(logFc) && abs(logFc) > log2FC
          inResult %>%
            mutate(valid = b)
        }) %>%
        .$valid

      finalResult <- if ((isInter && all(bs)) || (!isInter && any(bs))) {
        result
      }else {
        result %>%
          filter(FALSE)
      }
      finalResult
    }) %>%
    ungroup()
} else {
  finalDf
}

write.csv(finalDf, "finalDf.csv", row.names = FALSE)

finalDf <- if (execTrend) {
  finalDf %>%
    rowwise() %>%
    do({
      result <- as_tibble(.)
      method <- result$test.method %>%
        as.character()
      cMethod <- if (startsWith(method, "kruskal.test")) {
        "Median"
      }else "Mean"
      groups <- groupSort %>%
        str_split(",") %>%
        unlist()

      colNames <- groups %>%
        map_chr(function(group) {
          str_c(group, ".", cMethod)
        })
      data <- result %>%
        select(colNames) %>%
        unlist()
      increaseData <- sort(data)
      decreaseData <- sort(data, decreasing = T)
      b <- (trendMethod == "increase" & all(increaseData == data)) |
        (trendMethod == "decrease" & all(decreaseData == data)) |
        (trendMethod == "both" & (all(increaseData == data) | all(decreaseData == data)))
      finalResult <- if (b) {
        result
      }else {
        result %>%
          filter(FALSE)
      }
      finalResult
    }) %>%
    ungroup()
} else {
  finalDf
}

outFileName <- paste0(parent, "/Diff_Met_Test.csv")
write_csv(finalDf, outFileName)

diffData <- finalDf %>%
  select(c("Metabolite"))

write_csv(diffData, "Diff_Metabolite.csv")


